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Abstract. Recently, a whole class of evergy-preservingintegrators has been derived for the numerical solution of Hamiltonian 
problems 0,0,0]. In the mainstream of this research we have defined a new family of symplectic integrators depending 
on a real parameter a 0]. For a = 0, the corresponding method in the family becomes the classical Gauss collocation formula 
of order 2s, where s denotes the number of the internal stages. For any given non-null a, the corresponding method remains 
symplectic and has order 2s — 2: hence it may be interpreted as a 0(h 2s ~ 2 ) (symplectic) perturbation of the Gauss method. 
Under suitable assumptions, it can be shown that the parameter a may be properly tuned, at each step of the integration 
procedure, so as to guarantee energy conservation in the numerical solution. The resulting method shares the same order 2s 
as the generating Gauss formula, and is able to preserve both energy and quadratic invariants. 
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(/ is the identity matrix of dimension ni), two main lines of investigation may be traced in the current literature, having 
as objective the definition and the study of symplectic methods and energy-conserving methods, respectively. In fact, 
the symplecticity of the map and the conservation of the energy function are the most relevant features characterizing 
a Hamiltonian system. 

From the very beginning of this research activity, high order symplectic formulae were already available within 
the class of Runge-Kutta methods fill l24l |23ll . the Gauss collocation formulae being one noticeable example. One 
important implication of symplecticity of the discrete flow, for Gauss-Legendre methods, is the conservation of 
quadratic invariants. This circumstance makes the symplecticity property of the method particularly appealing in 
the numerical simulation of isolated mechanical systems in the form (Q~|i, since it provides a precise conservation 
of the total angular momentum during the time evolution of the state vector. As a further positive consequence, 
a symplectic method also conserves quadratic Hamiltonian functions (see the monographs IU5l l20ll for a thorough 
analysis of symplectic methods). 

Conversely, if one excludes the quadratic case, energy-conserving methods were initially not known within classical 
integration methods. The unsuccessful attempts to derive energy-preserving Runge-Kutta methods for polynomial 
Hamiltonians, culminated in the general feeling that such methods could not even exist (see lfl9ll and U). A completely 
new approach is represented by discrete gradient methods which are based upon the definition of a discrete counterpart 
of the gradient operator so that energy conservation of the numerical solution is guaranteed at each step and whatever 
the choice of the stepsize of integration (see lfl3ll2llo . 

More recently, the conservation of ener gy h as been approached by means of the definition of the discrete line in- 
tegral, in a series of papers (such as lfl7l Il8ll ). leading to the definition of Hamiltonian Boundary Value Methods 
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INTRODUCTION 



When dealing with the numerical integration of canonical Hamiltonian systems in the form 




(1) 



(HBVMs) (see for example HUHUQl]). They 

are a class of methods able to preserve, in the discrete solution, poly- 
nomial Hamiltonians of arbitrarily high degree (and, hence, a practical conservation of any sufficiently differentiable 
Hamiltonian. Such methods admit a Runge-Kutta formulation which reveals their close relationship with classical col- 
location formulae @]. An infinity extension of HBVMs has also been proposed in |5|] an d ITill . These limit methods 
may be interpreted as a generalization of the averaged vector field method defined in 12211 . 

Attempts to incorporate both symplecticity and energy conservation into the numerical method will clash with 
two non-existence results. The first fll2fl refers to non-integrable systems, that is systems that do not admit other 
independent first integrals different from the Hamiltonian function itself. According to the authors' words, it states 
that 

If [the method] is symplectic, and conserved H exactly, then it is the time advance map for the exact 
Hamiltonian system up to a reparametrization of time. 

The second negative result fioll refers to B-series symplectic methods applied to general (not necessarily non- 
integrable) Hamiltonian systems: 

The only symplectic method ( as B-series) that conserves the Hamiltonian for arbitrary H(y) is the exact flow 
of the differential equation. 

Despite these discouraging results, in |0] a new class of symplectic integrators of arbitrarily high-order has been 
proposed which, under some mild assumptions (see the next section), may share both features, in the sense specified 
in the theorem below. We prefer the use of the term "integrator" rather than method since, strictly speaking, our 
integrator may select a different symplectic formula from one integration step to the next, in order to enforce the energy 
conservation property. In what follows, we sketch the main ideas behind this approach. For further generalizations, 
as well as for a number of numerical evidences, we refer to ]7[]. We will begin with introducing a family of one-step 
methods 

yi(a,h) = & h (y ,a) (2) 
(h is the stepsize of integration), depending on a real parameter a, with the following specifics: 

1 . for any fixed choice of a ^ 0, the corresponding method is a symplectic Runge-Kutta method with s stages and 
of order 2s — 2, which exactly conserves all quadratic invariants; 

2. for a = one gets the Gauss collocation method (of order 2s); 

3. for any choice of yo and in a given range of the stepsize h, there exists a value of the parameter, say Oo, depending 
on yo and h, such that H{y\ (ao,h)) = H(yo) (energy conservation). 

The parametric method ([2]) realizes a symplectic perturbation of the Gauss method of size 0(h 2s ~ 2 ). Under suitable 
assumptions, as the parameter a ranges in a small interval centered at zero, the value of the numerical Hamiltonian 
function H(y\ ) will match H(y(to + h)) thus leading to energy conservation. This result is formalized as follows: 

Theorem 1 (Energy conservation) Under suitable assumptions, there exists a real sequence {cfy} such that the 
numerical solution defined by yk+i — &h{yki OCt), with yo defined in (Q, satisfies H(y/ C ) = H(yo). 

One important remark is in order to clarify this statement and how it relates to the above non-existence results. Let 
us select the value of the parameter a = Oq, if any, in order to enforce the energy conservation between the two state 
vectors yo and yi, as indicated at item 3 above@: the map y n> O/^y, (Xq) is symplectic and, by definition, assures the 
energy conservation condition H(y\ ) = H(yo). However, it is worth noticing that it would fail to provide a conservation 
of the Hamiltonian function if we changed the initial condition yo or the stepsize h. For example, in general for any 
So yo, we would obtain //(<!>/, (yo, <%)) ^ H(yo): in this case we should change the value of the parameter a in 
order to recover the equality condition^ Strictly speaking, the energy conservation property described in TheoremQ] 
weakens the standard energy conservation condition mentioned in the two non-existence results stated above and hence 
our methods are not meant to produce a counterexample of these statements. 



2 To avoid any misunderstanding, we emphasize that the value ceo is now maintained constant, otherwise the map would fail to be symplectic. 

3 More in general, the sequence {&/,} that will satisfy Theorem[T]starting at fo will differ from the sequence {oik}. Such sequences will be defined 
as the solution of the nonlinear system (7), as described in the next section. 



DEFINITION OF THE METHODS 



Let {c\ < Co < ■ ■ ■ < c v } and {b\, . . . ,b s } be the abscissae and the weights of the Gauss-Legendre quadrature 
formula in the interval [0,1]. We consider the Legendre polynomials Pj{t) of degree j — for j = l,...,s, shifted and 
normalized in the interval [0, 1] so that Jg Pi{x)Pj{x)dx = for i,j = 1, ... ,s, (5,y is the Kronecker symbol), and the 
sxs matrix 3 s = (P/(c,)). Our starting point is the following well-known decomposition of the Butcher array A of the 
Gauss method of order 2s M PP. 77-84]: 

A=3>X S 3>- 1 , 



where X s is defined as 
/ 



\ 



4- 







with 



4 = 



2\/47 



./' = l,...,s- 1. 



(3) 



(4) 



We now consider the matrix X s (a) obtained by perturbing (O as follows: 



X s {a) = 



( \ 
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-4i 
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-(4-i 




■a) 



= X s + aW s , 



(5) 



where a is a real parameter, W s = (e s e T sX — e s _\e T s \ and, as usual, 



methods (O we are interested in, is formally defined by the following tableau (see ©-(O): 



is the y'th unit vector. The family of 
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with 



a) = &>X s {a)&~ A =A + a£ g W s 3 g 



-i 



(6) 



Therefore A (0) =A and, moreover, the following result holds true [01. 

Theorem 2 For any fixed value of a, the Runge-Kutta method is symmetric and symplectic. For a — 0, the usual 
Gauss-Legendre method of order 2s is recovered. For any fixed CC ^ 0, a method of order 2s — 2 is obtained. 

If we can choose a = Oq so that the energy-conservation property be satisfied, then, one obtains a (symmetric), 
Energy and Quadratic Invariants Preserving (EQUIP) method, as specified in TheoremQ] of Gaussian type. Indeed, 
the conservation of quadratic invariants easily follows from the structure of the matrix © defining the method. In 
conclusion, these methods will provide an exact conservation of all quadratic invariants, besides the Hamiltonian 
function. In more details, if we denote, as usual, Y = (Y[ . . . Y[) T the vector of the stages, e = (1, . . . , l) T £ W, and 
defining the error function g(a,h) = H(y\ (cc,h)) — H(yo), the nonlinear system, in the unknowns Y\,... ,Y S and a, 
that is to be solved at each step for getting energy conservation, reads, for the given stepsize h, 



(7) 



Y = e®y + h(jtf(a)®I)F(Y), 
g(a,h) = 0. 

Concerning the question about the existence of a solution of (UJi, we make the following assumptions: 
(jz/\) the function g is analytical in a rectangle [—a, a] x [—A, A] centered at the origin; 

(s^2) let d be the order of the error in the Hamiltonian function associated with the Gauss method applied to the given 
Hamiltonian system (HJ and the given state vector yo, that is: 

g(0,h) = H(yi(0,A)) -H(y ) = c h d + 0(h d+1 ), c Q ^ 0. 

Then, we assume that for, any fixed a ^ in a suitable neighborhood of the origin, 

g(a,h) =c{a)h d ~ 2 + 0{h. d - 1 ), c(a)^0. 



Remark 1 We observe that, excluding the case where the Hamiltonian H(y) is quadratic (which would imply 
g(a,h) =0, for all a ), the error in the numerical Hamiltonian function associated with the Gauss method is expected 
to behave as 0(h 2s+l ). Consequently, d > 2s. 

The following result then holds true J3l ■ 

Theorem 3 Under the assumptions ( si\) and ( si%), there exists a function ccq = Go(h), defined in a neighborhood of 
the origin (—ho, ho), such that: 

(i) g(do(h),h) = 0, for all h G (-ho, ho); (ii) <Xo(h) — const • h 2 + 0(h 3 ). 

The next result concerns the order of convergence of the method © (again, the proof can be found in |01). 

Theorem 4 Consider the parametric method (O and suppose that the parameter a is actually a function of the stepsize 
h, according to what stated in Theorem\3\ Then, the resulting method has order 2s. 

Numerical tests concerning the new EQUIP methods of Gaussian type can be found in [7] and in the companion 
paper @]. 
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